% Copyright (C) 2019-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer
%
% This is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% It is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% For a copy of the GNU General Public License,
% see <http://www.gnu.org/licenses/>.

%% Generates:
% - Table B.2: First and second moments in the two exchange-rate regimes
% - Figure B.1: Quantitative model – ergodic debt distribution

rng('default')

% Display moments of float and peg
addpath('../Auxiliary_functions/')
clear;
if ~isfolder('Latex')
    mkdir('.','Latex');
end

labels_string=char('Unemployment rate','Consumption','Non-tr. cons.','Non-tr. output.','Trade balance','Real wage','Traded output','Interest rate','External debt','Debt-to-output ratio (ann)','Government spending share');
labels_string_tex=char('\bar{h}-h_t','c_t','c_t^N','y_t^N','y^T_t-c^T_t','w_t','y^T_t','r_t^{ann}','d_t','d_t/4(y^T_t+p_t^Nc^N_t)','G/Y');

T_ergodic = 1e6; %periods for simulation of mean;
T_burnin = 1e+3; %burn-in period

%load policy function
loaded_PFI_peg=load('Policy_functions/pfi_peg_g_low_beta_large_grid.mat');

if ~isfield(loaded_PFI_peg.par,'phi_e')
    loaded_PFI_peg.par.phi_e= 0; %exchange rate policy parameter
end


%% Set initial conditions of exogenous states, debt, and past wage equal to their respective unconditional means
temp=unique(loaded_PFI_peg.grid.y_level);
y_0 = temp(round(length(temp)/2));
temp=unique(loaded_PFI_peg.grid.r_level);
r_0 = temp(round(length(temp)/2));
temp=unique(loaded_PFI_peg.grid.g_level);
g_0 = temp(round(length(temp)/2));
d_0=loaded_PFI_peg.grid.d_grid(round(length(loaded_PFI_peg.grid.d_grid)/2));
p_0=2;
w_0=loaded_PFI_peg.grid.w_grid_vec(round(length(loaded_PFI_peg.grid.w_grid_vec)/2));

%P_trans is the transition probability matrix of tradable output. Cumul_P_trans is the cumulative probability matrix (useful for drawing tradable output realizations)
Cum_prob = cumsum(loaded_PFI_peg.MC.P_trans,2);

[Y_base_peg, var_name, pos]=simulate_model_general(y_0,r_0,g_0,d_0,w_0,p_0,T_ergodic,loaded_PFI_peg.dp,loaded_PFI_peg.h,loaded_PFI_peg.grid,loaded_PFI_peg.par,Cum_prob,1);
Y_base_peg=Y_base_peg(T_burnin+1:end,:);
Y_base_peg_mean=mean(Y_base_peg);

simulation_moments_peg=[mean(Y_base_peg(:,pos.unemployment))      std(Y_base_peg(:,pos.unemployment)) 
                mean(Y_base_peg(:,pos.c))                   std(Y_base_peg(:,pos.c)./Y_base_peg_mean(:,pos.c))
                mean(Y_base_peg(:,pos.c_N))                 std(Y_base_peg(:,pos.c_N)./Y_base_peg_mean(:,pos.c_N))
                mean(Y_base_peg(:,pos.y_N))                 std(Y_base_peg(:,pos.y_N)./Y_base_peg_mean(:,pos.y_N))
                mean(Y_base_peg(:,pos.trade_balance))       std(Y_base_peg(:,pos.trade_balance))
                mean(Y_base_peg(:,pos.w))                   std(Y_base_peg(:,pos.w)./Y_base_peg_mean(:,pos.w))
                mean(Y_base_peg(:,pos.y_T))                 std(Y_base_peg(:,pos.y_T)./Y_base_peg_mean(:,pos.y_T))
                mean(Y_base_peg(:,pos.r_ann))               std(Y_base_peg(:,pos.r_ann))
                mean(Y_base_peg(:,pos.D))                   std(Y_base_peg(:,pos.D)./Y_base_peg_mean(:,pos.D))
                mean(Y_base_peg(:,pos.D_y_annual))          std(Y_base_peg(:,pos.D_y_annual))                
                mean(Y_base_peg(:,pos.g_share))             std(Y_base_peg(:,pos.g_share))                
                ];

debt_density_peg=NaN(loaded_PFI_peg.grid.n_d,1);
d_grid_size=loaded_PFI_peg.grid.d_grid(2)-loaded_PFI_peg.grid.d_grid(1);
for grid_iter=1:loaded_PFI_peg.grid.n_d
    debt_density_peg(grid_iter,1) = mean(Y_base_peg(T_burnin+1:end,pos.D)==loaded_PFI_peg.grid.d_grid(grid_iter))./d_grid_size; %prob. distribution of debt
end

%% LaTex table
headers_string=char(' ','Mean','Std');
dyntable('moments_peg',headers_string,labels_string,simulation_moments_peg,size(labels_string,2)+2,5,3)

dyn_latex_table_modified('moments_peg','Latex/moments_peg',headers_string,labels_string_tex,simulation_moments_peg,size(labels_string,2)+2,5,3);

% clear('loaded_PFI_peg');

%% float
loaded_PFI_float=load('Policy_functions/vfi_float_g_low_beta_large_grid.mat');
Cum_prob = cumsum(loaded_PFI_float.MC.P_trans,2);

[Y_base_float, var_name, pos]=simulate_model_float(y_0,r_0,g_0,d_0,w_0,p_0,T_ergodic,loaded_PFI_float.dp,loaded_PFI_float.grid,loaded_PFI_float.par,Cum_prob,1);
Y_base_float=Y_base_float(T_burnin+1:end,:);
Y_base_float_mean=mean(Y_base_float);

simulation_moments_float=[mean(Y_base_float(:,pos.unemployment))      std(Y_base_float(:,pos.unemployment)) 
                mean(Y_base_float(:,pos.c))                   std(Y_base_float(:,pos.c)./Y_base_float_mean(:,pos.c))
                mean(Y_base_float(:,pos.c_N))                 std(Y_base_float(:,pos.c_N)./Y_base_float_mean(:,pos.c_N))
                mean(Y_base_float(:,pos.y_N))                 std(Y_base_float(:,pos.y_N)./Y_base_float_mean(:,pos.y_N))
                mean(Y_base_float(:,pos.trade_balance))       std(Y_base_float(:,pos.trade_balance))
                mean(Y_base_float(:,pos.w))                   std(Y_base_float(:,pos.w)./Y_base_float_mean(:,pos.w))
                mean(Y_base_float(:,pos.y_T))                 std(Y_base_float(:,pos.y_T)./Y_base_float_mean(:,pos.y_T))
                mean(Y_base_float(:,pos.r_ann))               std(Y_base_float(:,pos.r_ann))
                mean(Y_base_float(:,pos.D))                   std(Y_base_float(:,pos.D)./Y_base_float_mean(:,pos.D))
                mean(Y_base_float(:,pos.D_y_annual))          std(Y_base_float(:,pos.D_y_annual))                
                mean(Y_base_float(:,pos.g_share))             std(Y_base_float(:,pos.g_share))                
                ];

debt_density_float=NaN(loaded_PFI_float.grid.n_d,1);
d_grid_size=loaded_PFI_float.grid.d_grid(2)-loaded_PFI_float.grid.d_grid(1);
for grid_iter=1:loaded_PFI_float.grid.n_d
    debt_density_float(grid_iter,1) = mean(Y_base_float(T_burnin+1:end,pos.D)==loaded_PFI_float.grid.d_grid(grid_iter))./d_grid_size; %prob. distribution of debt
end

%% LaTex table
headers_string=char(' ','Mean','Std');
dyntable('moments_float',headers_string,labels_string,simulation_moments_float,size(labels_string,2)+2,5,3)

dyn_latex_table_modified('moments','Latex/moments_float',headers_string,labels_string_tex,simulation_moments_float,size(labels_string,2)+2,5,3);


%% output
            
hh=figure('Name','Debt Distribution');
plot(loaded_PFI_peg.grid.d_grid,debt_density_peg,'-',loaded_PFI_float.grid.d_grid,debt_density_float,':','linewidth',4);

xlabel('External Debt (level)')
ylabel('Density')

xlim([loaded_PFI_peg.grid.d_grid(1) loaded_PFI_peg.grid.d_grid(end)]);
ylim([0 max([debt_density_peg;debt_density_float])]);
legend('Peg','Float','Location','NorthWest')
print('Figures/Debt_distribution','-depsc2')
saveas(hh,'Figures/Debt_distribution')


%% LaTex table
headers_string=char(' ','Mean (peg)','Std (peg)','Mean (float)','Std (float)');
dyntable('moments',headers_string,labels_string,[simulation_moments_peg simulation_moments_float],size(labels_string,2)+2,5,3)

dyn_latex_table_modified('moments','Latex/moments',headers_string,labels_string_tex,[simulation_moments_peg simulation_moments_float],size(labels_string,2)+2,5,3);